1 Introduction

The Rise 2 Flow study used a number of instruments during the entry, exit, and weekly sessions (additionally personal data was collected during the entry session).

For entry and exit sessions six instruments were used:

  • Mindful attention awareness scale (MAAS).1
  • Scale of positive and negative experience (SPANE).2
  • Personal well-being index (PWB).3
  • Positive thinking (PST).4
  • Self efficacy (SE).5
  • Perceived productivity (PP).6

In addition, a seventh instrument, mini-IPIP, was added to the exit session to control for personality type.7

For the weekly sessions the WHO-5 well-being index was used (WHO-5).8

To summarize, the study was executed according to,

  1. An entry survey, at time \(t_0\) (MAAS, etc.)
  2. A weekly survey, at time \(t_{{t_0 + 1}}, \ldots, t_{{t_0+13}}\)
  3. Daily surveys
  4. An exit survey, at time \(t_{0+11}\) (MAAS, etc. and IPIP)

The main question the study tries to answer is: Is there a difference in the responses (i.e., higher/lower values on responses), across all or among some instruments, over time?

In order to answer that question we can first of all compare outcomes of \(t_0\) and \(t_1\), and contrast with personal data, e.g., age, gender, IPIP. Next, the weekly session can be used to investigate the trend over time, and perhaps see where a noticeable effect, if it exists, shows up.

2 Data cleaning

First we load all five datasets. For each dataset, make sure we set correct variable types.

entry <- read.xlsx(here("data/R2F Quant Comp Data Entry.xlsx")) # all instruments
exit <- read.xlsx(here("data/R2F Quant Comp Data Exit.xlsx")) # same here
weekly <- read.xlsx(here("data/R2F Quant Weekly.xlsx"), detectDates = T)
ipip <- read.xlsx(here("data/R2F IPIP.xlsx"))
pers <- read.xlsx(here("data/R2F Pers Data.xlsx")) # personal data
daily <- read.xlsx(here("data/R2F Daily.xlsx"), detectDates = T)

2.1 entry and exit datasets

We need to set a number of columns as factors (ordered factors in many cases). Additionally, when we have Yes/No answers, as we do in some questions below, we set them to 0/1 since they will be modeled with a Bernoulli likelihood.

First we clean the entry dataset.

for(j in 3:17){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Almost never", "Very infrequently", "Somewhat infrequently", "Somewhat frequently", "Very frequently", "Almost always"), ordered = TRUE))
}

for(j in 18:29){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Very rarely or never", "Rarely", "Sometimes", "Often", "Very often or always"), ordered = TRUE))
}

for(j in 30:37){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Strongly disagree", "Disagree", "Slightly disagree", "Mixed or neither agree nor disagree", "Slightly agree", "Agree", "Strongly agree"), ordered = TRUE))
}

for(i in 38:59)
  entry[,eval(i)] <- ifelse(entry[eval(i)] == 'Yes', 1, 0)

for(j in 60:69){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Hardly true","Rather true","Exactly true"), ordered = TRUE))
}

for(j in 70:76){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("None of the time","A little of the time","Some of the time","Most of the time","All of the time"), ordered = TRUE))
}

for(j in 77:79){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("1","2","3","4","5","6","7","8","9","10"), ordered = TRUE))
}

for(j in 80:80){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("You were a lot worse than other workers","You were a little worse than other workers","You were somewhat worse than other workers","You were about average","You were somewhat better than other workers","You were a little better than other workers","You were a lot better than other workers"), ordered = TRUE))
}

Next, we clean the exit dataset, which consists of the same variables as above.

for(j in 3:17){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Almost never", "Very infrequently", "Somewhat infrequently", "Somewhat frequently", "Very frequently", "Almost always"), ordered = TRUE))
}

for(j in 18:29){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Very rarely or never", "Rarely", "Sometimes", "Often", "Very often or always"), ordered = TRUE))
}

for(j in 30:37){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Strongly disagree", "Disagree", "Slightly disagree", "Mixed or neither agree nor disagree", "Slightly agree", "Agree", "Strongly agree"), ordered = TRUE))
}

for(i in 38:59)
  exit[,eval(i)] <- ifelse(exit[eval(i)] == 'Yes', 1, 0)


for(j in 60:69){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Hardly true","Rather true","Exactly true"), ordered = TRUE))
}

for(j in 70:76){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("None of the time","A little of the time","Some of the time","Most of the time","All of the time"), ordered = TRUE))
}

for(j in 77:79){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("1","2","3","4","5","6","7","8","9","10"), ordered = TRUE))
}

for(j in 80:80){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("You were a lot worse than other workers","You were a little worse than other workers","You were somewhat worse than other workers","You were about average","You were somewhat better than other workers","You were a little better than other workers","You were a lot better than other workers"), ordered = TRUE))
}

2.2 weekly dataset

The weekly dataset contains five columns of the same type, which we set as ordered categorical.

for(j in 4:8){
  set(weekly, i=NULL, j=j, value=factor(weekly[[j]], levels = c("At no time", "Less than half of the time", "Some of the time", "More than half of the time",  "Most of the time", "All of the time"), ordered = TRUE))
}

2.3 ipip dataset

The mini-IPIP instrument uses only five levels for all 20 questions (ordered categorical).

for(j in 3:22){
  set(ipip, i=NULL, j=j, value=factor(ipip[[j]], levels = c("Disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Agree"), ordered = TRUE))
}

2.4 pers dataset

In the pers dataset we have a mix of variable types we need to deal with. Age should be scaled, Gender and Occupation set to 1/2, and Living conditions set to factors (unordered).

# convert Age to numeric and standardize (note we have two NAs, so we suppress the warning we get)
suppressWarnings(pers$Age <- scale(as.numeric(pers$Age)))

# set 1/2 as indicators for gender
pers$Gender_0_1 <- ifelse(pers$Gender_0_1 == "Man/Transman", 1, 2)

# many categories exist for occupation, but let's simplify this so we can 
# use this variable for something. Indicate if they are students or not.
pers$Occupation <- ifelse(pers$Occupation == "Student", 1, 2)

# make living condition a categorical factor (unordered)
pers$Living_1_4 <- factor(pers$Living_1_4, levels = c("I live by myself", "I live in shared housing", "I live with a partner", "I live with my family"))

# remove the last three columns since we have >64% NAs
pers[,7:9] <- NULL

2.5 Missingness analysis

Check how many cells are NA.

table(is.na(entry))
## 
## FALSE  TRUE 
##  6862    98
table(is.na(exit))
## 
## FALSE  TRUE 
##  2684    36
table(is.na(weekly))
## 
## FALSE  TRUE 
##  2390     2
table(is.na(pers))
## 
## FALSE  TRUE 
##   519     3

Clearly we have missing data here, but it’s only a few percentages (\(0.6\)\(1.8\)%). Let’s use complete case analysis for now (i.e., remove all rows containing NAs). Later, if needed, we can model the missingness in a principled Bayesian way.

entry <- entry[complete.cases(entry), ]
exit <- exit[complete.cases(exit), ]
weekly <- weekly[complete.cases(weekly), ]
pers <- pers[complete.cases(pers), ]

Remember, the main question the study tries to answer is: Is there a difference in the responses (i.e., higher/lower values on responses), across all or among some instruments, over time? This means that for now we can discard of the weekly data and instead focus on the entry and exit instruments, together with personal data (i.e., age, gender, occupation, and living conditions).

In order for us to conduct inferences along this line of thinking, we will make use of dummy variable regression estimators (DVRE), which is numerically identical to deviation score estimators. The DVRE approach dummy encodes our time variable \(t\) and sets an index \(0/1\), where \(t_0 = 0\) and \(t_1 = 1\). In short, each subject (ID) will have two rows where one row is the entry instrument \(t = 0\), and one row has the exit instrument \(t = 1\). We’ll then see if our predictors are useful for predicting the outcome (i.e., is there a significant difference in the \(\beta\) estimators) and then also investigate if the outcome at \(t_0\) equals \(t_1\) (and if they aren’t the same, how much do they differ?)

# create new column t, and set to 0 or 1
entry$t <- 0
exit$t <- 1

# Remove all suffixes in both files so we have the same column names
colnames(entry) <- gsub('_en', '', colnames(entry))
colnames(exit) <- gsub('_ex', '', colnames(exit))

# also remove all dashes and underscores in column names
colnames(entry) <- gsub('_', '', colnames(entry))
colnames(exit) <- gsub('_', '', colnames(exit))
colnames(entry) <- gsub('-', '', colnames(entry))
colnames(exit) <- gsub('-', '', colnames(exit))
colnames(ipip) <- gsub('_', '', colnames(ipip))
colnames(pers) <- gsub('_', '', colnames(pers))
colnames(ipip) <- gsub('-', '', colnames(ipip))
colnames(pers) <- gsub('-', '', colnames(pers))

# add rows of entry and exit
d <- rbind(entry, exit)

# keep only rows if they show up twice, so we have a within-measurement of each 
# individual at time t_0 and t_1
d <- ddply(d, .(ID), function(x){
    if(nrow(x)==1){
        return(NULL)
    }
    else{
        return(x)
    }
})

table(d$t)
## 
##  0  1 
## 13 13

We have \(13\) individuals with entry and exit surveys completely filled out. Let’s also add the ipip and pers datasets, and then use only individuals with complete data.

d <- inner_join(d, pers, by = "ID")
d <- inner_join(d, ipip, by = "ID")

After joining these two datasets we have 12 subjects in the dataset and the total sample size is \(n=\) 24.

3 Model designs

3.1 Assumptions concerning the data-generation process

Next, in order to design an appropriate model, one of the things we need to make assumptions about is the underlying data-generative model, i.e., what model generated the type of data we have at our disposal. In short, there are four candidates: Cumulative, Continuation ratio, Stopping ratio, and Adjacent category. Above we opted for a default Cumulative. Let’s instead use LOO to create identical models and then use the same data for all models to compare them from an information theoretical perspective.

bf1 <- bf(mvbind(MAASQ116,
                 MAASQ216,
                 MAASQ316,
                 MAASQ416,
                 MAASQ516,
                 MAASQ616,
                 MAASQ716,
                 MAASQ816,
                 MAASQ916,
                 MAASQ1016,
                 MAASQ1116,
                 MAASQ1216,
                 MAASQ1316,
                 MAASQ1416,
                 MAASQ1516) ~ 
            1 )

m0 <- brm(bf1,
    family = cumulative,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 17.5 seconds.
## Chain 3 finished in 17.5 seconds.
## Chain 1 finished in 17.8 seconds.
## Chain 4 finished in 18.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 17.7 seconds.
## Total execution time: 18.5 seconds.
mcr <- brm(bf1,
    family = cratio,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 4.1 seconds.
## Chain 1 finished in 4.1 seconds.
## Chain 3 finished in 5.3 seconds.
## Chain 4 finished in 5.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 4.7 seconds.
## Total execution time: 5.6 seconds.
msr <- brm(bf1,
    family = sratio,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 5.2 seconds.
## Chain 4 finished in 6.2 seconds.
## Chain 1 finished in 6.3 seconds.
## Chain 3 finished in 6.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 6.1 seconds.
## Total execution time: 6.8 seconds.
mac <- brm(bf1,
    family = acat,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 8.5 seconds.
## Chain 2 finished in 8.5 seconds.
## Chain 4 finished in 8.6 seconds.
## Chain 3 finished in 9.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 8.6 seconds.
## Total execution time: 9.3 seconds.
m0 <- add_criterion(m0, criterion = "loo")
msr <- add_criterion(msr, criterion = "loo")
mcr <- add_criterion(mcr, criterion = "loo")
mac <- add_criterion(mac, criterion = "loo")

loo_compare(m0, mcr, msr, mac)
##     elpd_diff se_diff
## m0   0.0       0.0   
## mcr -2.7       3.1   
## msr -2.8       3.1   
## mac -3.1       1.6

Evidently, since \(\Delta\)SE is fairly large, in comparison to the relative difference in expected log pointwise predictive density (\(\Delta\)elpd), we can safely assume that there’s no significant difference between the likelihoods and, thus, opt for the standard approach, i.e., the Cumulative distribution, in this case represented by \(\mathcal{M}_0\).

3.2 weekly data

3.2.1 Prior and posterior predictive checks

We want to model the weekly instrument and how it varies over time \(t\). In this case, \(t\) will be modeled using a Gaussian Process and then each individual in the dataset will be modeled using a varying intercept.

weekly <- left_join(weekly, pers, by = "ID")

# make sure ID is a factor
weekly$ID <- as.factor(weekly$ID)

# Need to format data so that each person's measurements goes 
# from 1 ... n, according to order, then we use those numbers as temporal var
weekly <- weekly[order(weekly$ID),]
weekly$Seq <- sequence(tabulate(weekly$ID))
# Pick first question. We could model this as multivariate but it won't give us 
# anything extra
p <- get_prior(WQ1_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = weekly,
               family = cumulative())

p$prior[1] <- "normal(0,2)"
p$prior[8] <- "normal(0,5)"
p$prior[13] <- ""
p$prior[14] <- "inv_gamma(1.6, 0.1)"
p$prior[15] <- "weibull(2,1)"

WQ1gp <- brm(WQ1_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             sample_prior = "only",
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 chains, at most 8 in parallel...
## 
## Chain 1 finished in 0.2 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
## Chain 2 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.7 seconds.
pp_check(WQ1gp, type="bars", nsamples = 100)

The priors are fairly uniform on the outcome space. Let’s sample with data now.

WQ1gp <- brm(WQ1_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.9999, max_treedepth = 13),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 90.4 seconds.
## Chain 1 finished in 91.5 seconds.
## Chain 3 finished in 95.1 seconds.
## Chain 4 finished in 98.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 93.9 seconds.
## Total execution time: 99.0 seconds.
WQ2gp <- brm(WQ2_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.99, max_treedepth = 12),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 71.1 seconds.
## Chain 3 finished in 88.0 seconds.
## Chain 2 finished in 104.7 seconds.
## Chain 4 finished in 106.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 92.6 seconds.
## Total execution time: 106.8 seconds.
WQ3gp <- brm(WQ3_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.99),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 47.7 seconds.
## Chain 1 finished in 50.9 seconds.
## Chain 2 finished in 51.2 seconds.
## Chain 4 finished in 53.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 50.7 seconds.
## Total execution time: 53.3 seconds.
WQ4gp <- brm(WQ4_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.999),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 chains, at most 8 in parallel...
## 
## Chain 4 finished in 67.8 seconds.
## Chain 3 finished in 74.0 seconds.
## Chain 2 finished in 85.3 seconds.
## Chain 1 finished in 87.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 78.5 seconds.
## Total execution time: 89.0 seconds.
WQ5gp <- brm(WQ5_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID),
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.99, max_treedepth = 12),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 126.8 seconds.
## Chain 4 finished in 130.7 seconds.
## Chain 1 finished in 132.6 seconds.
## Chain 2 finished in 133.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 130.9 seconds.
## Total execution time: 134.0 seconds.

3.2.2 Diagnostics

models <- list(WQ1gp, WQ2gp, WQ3gp, WQ4gp, WQ5gp)

# Check divergences, tree depth, energy for each model
for(m in models) {
  rstan::check_hmc_diagnostics(eval(m)$fit)
}
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# Check rhat and neff
for (m in models) {
  if (max(bayesplot::rhat(eval(m)), na.rm = T) >= 1.01)
    print("Warning: Rhat >= 1.01")
  
  if (min(bayesplot::neff_ratio(eval(m)), na.rm = T) <= 0.2)
    print("Warning: ESS <= 0.2")
}

3.2.3 Parameter estimates

In short, very little happening over time according to the weekly instrument.

3.3 daily data

3.3.1 Prior and posterior predictive checks

Let’s redo this but for the daily dataset.

daily <- left_join(daily, pers, by = "ID")
daily$ID <- as.factor(daily$ID)
daily$date_s <- as.numeric(daily$date)
daily$Living14 <- as.numeric(daily$Living14)
daily <- daily[order(weekly$ID),]
daily$Seq <- sequence(tabulate(daily$ID))

p <- get_prior(resp ~ 1 + gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = daily, family = cumulative)
p$prior[1] <- "normal(0,2)"
p$prior[6] <- "normal(0,5)"
p$prior[17] <- "inv_gamma(1.6, 0.1)"
p$prior[18] <- "weibull(2,1)"

m_daily <- brm(resp ~ 1 + gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = daily, 
               family = cumulative,
               prior = p,
               sample_prior = "only",
               silent = TRUE, 
               refresh = 0)#, threads = threading(4))
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 0.3 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.9 seconds.
pp_check(m_daily, type = "bars", nsamples = 100)

m_daily <- brm(resp ~ 1 + gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = daily, 
               family = cumulative,
               prior = p,
               control = list(adapt_delta=0.99),
               silent = TRUE, 
               refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 67.4 seconds.
## Chain 4 finished in 68.3 seconds.
## Chain 3 finished in 96.6 seconds.
## Chain 1 finished in 103.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 83.8 seconds.
## Total execution time: 103.3 seconds.
pp_check(m_daily, type = "bars", nsamples = 100)

3.3.2 Diagnostics

rstan::check_hmc_diagnostics(m_daily$fit)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# check rhat and ESS
if(max(bayesplot::rhat(eval(m_daily)), na.rm=T) >= 1.01) 
  print("Warning: Rhat >=1.01")

if(min(bayesplot::neff_ratio(eval(m_daily)), na.rm=T) <= 0.2) 
  print("Warning: ESS <=0.2")

3.3.3 Parameter estimates

Let’s plot an interesting effect, which albeit is not significant has a tendency that’s interesting.

plot(conditional_effects(m_daily), plot = FALSE)[[2]] +
  xlim(1,2.5) +
  scale_x_continuous(name=NULL, 
                     breaks = c(1,2), 
                     label = c("Man/\n transman","Woman/\n transwoman"), 
                     expand = c(0.08, 0)) +
  scale_y_continuous(name=NULL, breaks = c(6,7)) + 
  theme_tufte()

Next, look at the trend.

plot(conditional_effects(m_daily), plot = FALSE)[[5]] +
  xlab("Day") + ylab("") +
  theme_tufte()

One things clear, the uncertainty at the end makes an inferences hard

3.4 Multivariate models with temporal variable

We want a model that uses personal data (e.g., age) as predictors. The idea is to later see if any of the predictors have predictive capacity for our six different survey instruments.

First, design a model with predictors from personal data and compare with our null model (\(\mathcal{M}_0\)), together with our indicator \(t\) for time, and a varying intercept ID that varies depending on subject.

bf1 <- bf(mvbind(MAASQ116,
                 MAASQ216,
                 MAASQ316,
                 MAASQ416,
                 MAASQ516,
                 MAASQ616,
                 MAASQ716,
                 MAASQ816,
                 MAASQ916,
                 MAASQ1016,
                 MAASQ1116,
                 MAASQ1216,
                 MAASQ1316,
                 MAASQ1416,
                 MAASQ1516) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1|c|ID))

m_pers <- brm(bf1,
    family = cumulative,
    data = d,
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 175.6 seconds.
## Chain 1 finished in 186.2 seconds.
## Chain 4 finished in 194.9 seconds.
## Chain 3 finished in 198.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 188.9 seconds.
## Total execution time: 199.1 seconds.
    #threads = threading(4))

m_pers <- add_criterion(m_pers, criterion = "loo")

# compare the new model with our null model, it should hopefully be better...
(l <- loo_compare(m0, m_pers))
##        elpd_diff se_diff
## m_pers   0.0       0.0  
## m0     -51.3      24.5

Calculating the credible interval (\(z_{95\%}=1.96\)) we can see that it doesn’t cross zero, i.e., adding these predictors makes the relative out of sample prediction capability jump up even though we add seven population-level predictors and one group-level predictor (Living14 contains three parameters to estimate).

l[2,1] + c(-1,1) * 1.96 * l[2,2]
## [1] -87.327351   2.811668

3.4.1 Prior and posterior predictive checks

Let’s now sample all models with the same predictors we used above (visual prior and posterior predictive checks were done for each model).

################################################################################
#
# MAAS
#
################################################################################
bf1 <- bf(mvbind(MAASQ116,
                 MAASQ216,
                 MAASQ316,
                 MAASQ416,
                 MAASQ516,
                 MAASQ616,
                 MAASQ716,
                 MAASQ816,
                 MAASQ916,
                 MAASQ1016,
                 MAASQ1116,
                 MAASQ1216,
                 MAASQ1316,
                 MAASQ1416,
                 MAASQ1516) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

# Prior predictive checks have been conducted and the probability mass 
# is distributed evenly on the outcome space
p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))
  
         
m_maas <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    #sample_prior = "only", # use if you only want to sample from prior
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 31.2 seconds.
## Chain 1 finished in 31.3 seconds.
## Chain 2 finished in 31.8 seconds.
## Chain 3 finished in 32.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 31.7 seconds.
## Total execution time: 32.6 seconds.
    #threads = threading(4))

# use the below to do prior and posterior predictive checks
# change resp to the response variable you want to check 
# pp_check(m_maas, resp = "MAASQ116", type = "bars, nsamples = 250) 
################################################################################
#
# SPANE
#
################################################################################

bf1 <- bf(mvbind(SPANEQ115,
                 SPANEQ215,
                 SPANEQ315,
                 SPANEQ415,
                 SPANEQ515,
                 SPANEQ615,
                 SPANEQ715,
                 SPANEQ815,
                 SPANEQ915,
                 SPANEQ1015,
                 SPANEQ1115,
                 SPANEQ1215) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_spane <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only", # use if you only want to sample from prior
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 21.5 seconds.
## Chain 4 finished in 21.7 seconds.
## Chain 2 finished in 21.9 seconds.
## Chain 1 finished in 22.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 21.8 seconds.
## Total execution time: 22.4 seconds.
    #threads = threading(4))

# use the below to do prior and posterior predictive checks
# change resp to the response variable you want to check 
# pp_check(m_spane, resp = "SPANEQ115", type = "bars", nsamples = 250) 
################################################################################
#
# PWB
#
################################################################################

bf1 <- bf(mvbind(PWBQ117,
                 PWBQ217,
                 PWBQ317,
                 PWBQ417,
                 PWBQ517,
                 PWBQ617,
                 PWBQ717,
                 PWBQ817) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_pwb <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only", # use if you only want to sample from prior
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 18.0 seconds.
## Chain 1 finished in 18.4 seconds.
## Chain 2 finished in 18.4 seconds.
## Chain 3 finished in 19.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 18.5 seconds.
## Total execution time: 19.5 seconds.
    #threads = threading(4))

# pp_check(m_pwb, resp = "PWBQ117", type = "bars", nsamples = 250) 
################################################################################
#
# PST
#
################################################################################
# Here we move to a bernoulli since we have 0/1 outcome
bf1 <- bf(mvbind(PSTQ101,
                 PSTQ201,
                 PSTQ301,
                 PSTQ401,
                 PSTQ501,
                 PSTQ601,
                 PSTQ701,
                 PSTQ801,
                 PSTQ901,
                 PSTQ1001,
                 PSTQ1101,
                 PSTQ1201,
                 PSTQ1301,
                 PSTQ1401,
                 PSTQ1501,
                 PSTQ1601,
                 PSTQ1701,
                 PSTQ1801,
                 PSTQ1901,
                 PSTQ2001,
                 PSTQ2101,
                 PSTQ2201) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=bernoulli) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_pst <- brm(bf1,
    family = bernoulli,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 11.0 seconds.
## Chain 3 finished in 11.1 seconds.
## Chain 2 finished in 11.2 seconds.
## Chain 4 finished in 11.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 11.2 seconds.
## Total execution time: 11.7 seconds.
# pp_check(m_pst, resp = "PSTQ101", type = "bars", nsamples = 250) 
################################################################################
#
# SE - two models (Cumulative and Bernoulli), both multivariate 
#
################################################################################
# NOTE: the first two questions can be analyzed as Bernoulli 
bf1 <- bf(mvbind(#SEQ114, # this and the next outcome only has two categories
                 #SEQ214,
                 SEQ314,
                 SEQ414,
                 SEQ514,
                 SEQ614,
                 SEQ714,
                 SEQ814,
                 SEQ914,
                 SEQ1014) ~
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_se <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 8.0 seconds.
## Chain 1 finished in 8.4 seconds.
## Chain 2 finished in 9.4 seconds.
## Chain 3 finished in 9.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 8.8 seconds.
## Total execution time: 9.8 seconds.
    #threads = threading(4))

# pp_check(m_seq, resp = "SEQ314", type = "bars", nsamples = 250) 

# next take the two questions that only had answers on two levels
bf1 <- bf(mvbind(SEQ114,
                 SEQ214) ~
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=bernoulli) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_se2 <- brm(bf1,
    family = bernoulli,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.7 seconds.
# pp_check(m_seq2, resp = "SEQ114", type = "bars", nsamples = 250) 
################################################################################
#
# PPHQ
#
################################################################################


bf1 <- bf(mvbind(PPHQ115,
                 PPHQ215,
                 PPHQ315,
                 PPHQ415,
                 PPHQ515,
                 PPHQ615,
                 PPHQ715,
                 PPRQ1110,
                 PPRQ2110,
                 PPRQ3110,
                 PPOQ117
                 ) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_pph <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0,
    # threads = threading(4)
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 18.5 seconds.
## Chain 2 finished in 18.8 seconds.
## Chain 4 finished in 19.2 seconds.
## Chain 1 finished in 19.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 19.0 seconds.
## Total execution time: 19.8 seconds.
# pp_check(m_pphq, resp = "PPHQ115", type = "bars", nsamples = 250) 

So we have now sampled seven models. For each instrument we have one multivariate model (\(>1\) outcome), where we also estimate the correlation between the outcome given the subject’s ID.

Most of the models are using, as we’ve already discussed, the Cumulative family. However, in a few cases we use the Bernoulli family. First, we use it for outcomes that were binary (i.e., Yes/No outcomes). Second, we use it for two questions in the SE instrument (\(\mathcal{M}_{\text{SE}2}\)), which only had empirical outcomes on two levels.

Before we put any trust in these models we need to check a number of diagnostics.

3.4.2 Diagnostics

First, we check general HMC diagnostics. Second, we check that \(\widehat{R} \leq 1.01\) and that the effective sample size (ESS) is \(\geq 0.2\).

models <- list(m_maas, m_spane, m_pwb, m_pst, m_se, m_se2, m_pph)

# Check divergences, tree depth, energy for each model
for(m in models) {
  rstan::check_hmc_diagnostics(eval(m)$fit)
}
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# Check rhat and neff
for (m in models) {
  if (max(bayesplot::rhat(eval(m)), na.rm = T) >= 1.01)
    # diag in corr matrix always NA
    print("Warning: Rhat >= 1.01")
  
  if (min(bayesplot::neff_ratio(eval(m)), na.rm = T) <= 0.2)
    print("Warning: ESS <= 0.2")
}

Hamiltonian Monte Carlo diagnostics indicate all is well (one can also check traceplots for all models, i.e., using plot(model)).

3.4.3 Parameter estimates

First, let’s summarize what we have so far:

  1. Weekly trends modeled with Gaussian Processes (results in Sect. 3.2).
  2. Daily trends modeled with Gaussian Processes (results in Sect. 3.3). We’ve also shown conditional effects of gender since it is the main driver for deviations.
  3. Seven multivariate models with temporal variable to analyze entry/exit. Additionally, we have an eight model for analyzing two questions (Q1 and Q2) in the SE instrument (using a Bernoulli likelihood).

We can now investigate which, if any, parameters are significant on the 95%-level for the third item in the list, i.e., for each model, and each question,plot all parameters. Any parameters that are significant could warrant further analysis.

For each model we first analyze the temporal variable \(t\). Then we analyze the rest of the \(\beta\) parameters. Finally, where appropriate we look at the conditional effects. All densities we plot are limited to 95% probability mass, while the 50% probability mass is shaded around a vertical line, which is the median.

Each of these models has \(n=\) 24.

3.4.3.1 Mindful attention awareness scale (MAAS)

The model estimated 1355 parameters.

For Questions \(1\)\(3\), \(8\), and \(14\) the parameter \(t\) is negative. In short, the respondents answered with lower responses at \(t_1\) than at \(t_0\) for these questions. If we look at each unique parameter, in each question, can we see which, if any, predictors drove these results?

Much variance makes the estimates uncertain and, ultimately, no parameter is significant (not even Q7, even though it looks that way).

Two notable parameters, first for Q1 (woman/transwoman answers significantly higher than man/transman), while in Q6 man/transman answers significantly lower.

Occupation status (student/non-student) does not have any effect. Let’s, finally, look at living conditions before we turn our attention to the next instrument.

For Q15, if living with a partner, the respondent will respond significantly higher. The same for Q13 if the respondent lives with their family.

3.4.3.2 Scale of positive and negative experience (SPANE)

3.4.3.3 Personal well-being index (PWB)

For Q4 the respondents gave higher responses at \(t_1\). Let’s look at the parameters for that question.

No parameters influencing that question to a larger extent. In short, respondent, no matter their age, etc., responded higher at \(t_1\).

3.4.3.4 Positive thinking (PST)

3.4.3.5 Self efficacy (SE)

3.4.3.6 Perceived productivity (PP)

Questions \(2\), \(6\), and \(11\), show a significant difference when moving from \(t_0\) to \(t_1\). Let’s examine the parameters for these three questions.

Age has a negative impact (lower age lower on the response), however it is not significant.

##  Family: MV(cumulative, cumulative, cumulative, cumulative, cumulative, cumulative, cumulative, cumulative, cumulative, cumulative, cumulative) 
##   Links: mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity
##          mu = logit; disc = identity 
## Formula: PPHQ115 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPHQ215 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPHQ315 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPHQ415 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPHQ515 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPHQ615 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPHQ715 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPRQ1110 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPRQ2110 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPRQ3110 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##          PPOQ117 ~ 1 + Age + Gender01 + Occupation + Living14 + t + (1 | c | ID) 
##    Data: d (Number of observations: 24) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~ID (Number of levels: 12) 
##                                            Estimate Est.Error l-95% CI u-95% CI
## sd(PPHQ115_Intercept)                          1.41      0.50     0.47     2.37
## sd(PPHQ215_Intercept)                          1.68      0.47     0.80     2.60
## sd(PPHQ315_Intercept)                          1.18      0.48     0.25     2.15
## sd(PPHQ415_Intercept)                          1.51      0.47     0.64     2.47
## sd(PPHQ515_Intercept)                          1.25      0.46     0.38     2.15
## sd(PPHQ615_Intercept)                          1.28      0.45     0.44     2.22
## sd(PPHQ715_Intercept)                          0.81      0.42     0.14     1.76
## sd(PPRQ1110_Intercept)                         1.00      0.45     0.21     1.93
## sd(PPRQ2110_Intercept)                         1.48      0.46     0.61     2.36
## sd(PPRQ3110_Intercept)                         1.14      0.43     0.37     2.04
## sd(PPOQ117_Intercept)                          1.04      0.46     0.21     1.95
## cor(PPHQ115_Intercept,PPHQ215_Intercept)      -0.23      0.23    -0.64     0.24
## cor(PPHQ115_Intercept,PPHQ315_Intercept)      -0.12      0.24    -0.58     0.36
## cor(PPHQ215_Intercept,PPHQ315_Intercept)       0.20      0.24    -0.30     0.64
## cor(PPHQ115_Intercept,PPHQ415_Intercept)      -0.12      0.24    -0.56     0.35
## cor(PPHQ215_Intercept,PPHQ415_Intercept)       0.24      0.23    -0.22     0.65
## cor(PPHQ315_Intercept,PPHQ415_Intercept)       0.17      0.25    -0.34     0.62
## cor(PPHQ115_Intercept,PPHQ515_Intercept)      -0.11      0.25    -0.59     0.37
## cor(PPHQ215_Intercept,PPHQ515_Intercept)       0.26      0.23    -0.22     0.66
## cor(PPHQ315_Intercept,PPHQ515_Intercept)       0.13      0.24    -0.35     0.60
## cor(PPHQ415_Intercept,PPHQ515_Intercept)       0.22      0.24    -0.29     0.65
## cor(PPHQ115_Intercept,PPHQ615_Intercept)      -0.23      0.24    -0.67     0.26
## cor(PPHQ215_Intercept,PPHQ615_Intercept)       0.26      0.22    -0.21     0.67
## cor(PPHQ315_Intercept,PPHQ615_Intercept)       0.17      0.26    -0.35     0.64
## cor(PPHQ415_Intercept,PPHQ615_Intercept)       0.22      0.24    -0.28     0.65
## cor(PPHQ515_Intercept,PPHQ615_Intercept)       0.18      0.24    -0.32     0.61
## cor(PPHQ115_Intercept,PPHQ715_Intercept)      -0.03      0.26    -0.51     0.46
## cor(PPHQ215_Intercept,PPHQ715_Intercept)       0.08      0.26    -0.44     0.56
## cor(PPHQ315_Intercept,PPHQ715_Intercept)       0.02      0.26    -0.51     0.52
## cor(PPHQ415_Intercept,PPHQ715_Intercept)       0.08      0.25    -0.43     0.56
## cor(PPHQ515_Intercept,PPHQ715_Intercept)       0.09      0.25    -0.40     0.56
## cor(PPHQ615_Intercept,PPHQ715_Intercept)       0.04      0.26    -0.46     0.55
## cor(PPHQ115_Intercept,PPRQ1110_Intercept)      0.06      0.25    -0.44     0.53
## cor(PPHQ215_Intercept,PPRQ1110_Intercept)     -0.10      0.24    -0.53     0.39
## cor(PPHQ315_Intercept,PPRQ1110_Intercept)     -0.03      0.26    -0.53     0.48
## cor(PPHQ415_Intercept,PPRQ1110_Intercept)      0.04      0.25    -0.44     0.51
## cor(PPHQ515_Intercept,PPRQ1110_Intercept)     -0.04      0.25    -0.51     0.45
## cor(PPHQ615_Intercept,PPRQ1110_Intercept)     -0.03      0.25    -0.52     0.45
## cor(PPHQ715_Intercept,PPRQ1110_Intercept)     -0.03      0.26    -0.52     0.46
## cor(PPHQ115_Intercept,PPRQ2110_Intercept)      0.17      0.24    -0.32     0.61
## cor(PPHQ215_Intercept,PPRQ2110_Intercept)     -0.32      0.22    -0.70     0.16
## cor(PPHQ315_Intercept,PPRQ2110_Intercept)     -0.10      0.24    -0.54     0.40
## cor(PPHQ415_Intercept,PPRQ2110_Intercept)     -0.19      0.23    -0.62     0.28
## cor(PPHQ515_Intercept,PPRQ2110_Intercept)     -0.24      0.24    -0.66     0.25
## cor(PPHQ615_Intercept,PPRQ2110_Intercept)     -0.19      0.23    -0.62     0.28
## cor(PPHQ715_Intercept,PPRQ2110_Intercept)     -0.11      0.25    -0.59     0.40
## cor(PPRQ1110_Intercept,PPRQ2110_Intercept)     0.10      0.25    -0.41     0.56
## cor(PPHQ115_Intercept,PPRQ3110_Intercept)      0.17      0.25    -0.34     0.61
## cor(PPHQ215_Intercept,PPRQ3110_Intercept)     -0.26      0.23    -0.66     0.22
## cor(PPHQ315_Intercept,PPRQ3110_Intercept)     -0.13      0.25    -0.59     0.36
## cor(PPHQ415_Intercept,PPRQ3110_Intercept)     -0.13      0.24    -0.58     0.35
## cor(PPHQ515_Intercept,PPRQ3110_Intercept)     -0.16      0.25    -0.61     0.34
## cor(PPHQ615_Intercept,PPRQ3110_Intercept)     -0.17      0.24    -0.63     0.32
## cor(PPHQ715_Intercept,PPRQ3110_Intercept)     -0.07      0.26    -0.54     0.46
## cor(PPRQ1110_Intercept,PPRQ3110_Intercept)     0.12      0.25    -0.38     0.58
## cor(PPRQ2110_Intercept,PPRQ3110_Intercept)     0.23      0.24    -0.27     0.66
## cor(PPHQ115_Intercept,PPOQ117_Intercept)      -0.01      0.26    -0.51     0.48
## cor(PPHQ215_Intercept,PPOQ117_Intercept)      -0.14      0.24    -0.59     0.35
## cor(PPHQ315_Intercept,PPOQ117_Intercept)      -0.11      0.25    -0.57     0.40
## cor(PPHQ415_Intercept,PPOQ117_Intercept)      -0.18      0.24    -0.64     0.32
## cor(PPHQ515_Intercept,PPOQ117_Intercept)      -0.15      0.24    -0.59     0.35
## cor(PPHQ615_Intercept,PPOQ117_Intercept)      -0.07      0.25    -0.54     0.41
## cor(PPHQ715_Intercept,PPOQ117_Intercept)      -0.06      0.26    -0.56     0.45
## cor(PPRQ1110_Intercept,PPOQ117_Intercept)     -0.02      0.26    -0.50     0.48
## cor(PPRQ2110_Intercept,PPOQ117_Intercept)      0.11      0.25    -0.38     0.57
## cor(PPRQ3110_Intercept,PPOQ117_Intercept)      0.08      0.26    -0.44     0.55
##                                            Rhat Bulk_ESS Tail_ESS
## sd(PPHQ115_Intercept)                      1.00     1780     1119
## sd(PPHQ215_Intercept)                      1.00     2612     1982
## sd(PPHQ315_Intercept)                      1.00     1863     1394
## sd(PPHQ415_Intercept)                      1.00     2661     2285
## sd(PPHQ515_Intercept)                      1.00     2330     1648
## sd(PPHQ615_Intercept)                      1.00     2900     1950
## sd(PPHQ715_Intercept)                      1.00     2362     2261
## sd(PPRQ1110_Intercept)                     1.00     2162     1951
## sd(PPRQ2110_Intercept)                     1.00     2484     1634
## sd(PPRQ3110_Intercept)                     1.00     1946     1414
## sd(PPOQ117_Intercept)                      1.00     2142     1377
## cor(PPHQ115_Intercept,PPHQ215_Intercept)   1.00     2411     2234
## cor(PPHQ115_Intercept,PPHQ315_Intercept)   1.00     4034     2772
## cor(PPHQ215_Intercept,PPHQ315_Intercept)   1.00     3428     3035
## cor(PPHQ115_Intercept,PPHQ415_Intercept)   1.00     3267     3303
## cor(PPHQ215_Intercept,PPHQ415_Intercept)   1.00     3727     3342
## cor(PPHQ315_Intercept,PPHQ415_Intercept)   1.00     3099     2770
## cor(PPHQ115_Intercept,PPHQ515_Intercept)   1.00     3213     2828
## cor(PPHQ215_Intercept,PPHQ515_Intercept)   1.00     3376     2844
## cor(PPHQ315_Intercept,PPHQ515_Intercept)   1.00     3746     3316
## cor(PPHQ415_Intercept,PPHQ515_Intercept)   1.00     3510     3397
## cor(PPHQ115_Intercept,PPHQ615_Intercept)   1.00     3393     2787
## cor(PPHQ215_Intercept,PPHQ615_Intercept)   1.00     3648     3067
## cor(PPHQ315_Intercept,PPHQ615_Intercept)   1.00     4039     3109
## cor(PPHQ415_Intercept,PPHQ615_Intercept)   1.00     3675     2785
## cor(PPHQ515_Intercept,PPHQ615_Intercept)   1.00     3738     3278
## cor(PPHQ115_Intercept,PPHQ715_Intercept)   1.00     5558     3167
## cor(PPHQ215_Intercept,PPHQ715_Intercept)   1.00     4826     3271
## cor(PPHQ315_Intercept,PPHQ715_Intercept)   1.00     4665     3041
## cor(PPHQ415_Intercept,PPHQ715_Intercept)   1.00     4318     3318
## cor(PPHQ515_Intercept,PPHQ715_Intercept)   1.00     4054     3026
## cor(PPHQ615_Intercept,PPHQ715_Intercept)   1.00     3974     3156
## cor(PPHQ115_Intercept,PPRQ1110_Intercept)  1.00     4149     3586
## cor(PPHQ215_Intercept,PPRQ1110_Intercept)  1.00     4025     3102
## cor(PPHQ315_Intercept,PPRQ1110_Intercept)  1.00     4067     3243
## cor(PPHQ415_Intercept,PPRQ1110_Intercept)  1.00     4133     3153
## cor(PPHQ515_Intercept,PPRQ1110_Intercept)  1.00     4135     3150
## cor(PPHQ615_Intercept,PPRQ1110_Intercept)  1.00     3402     3392
## cor(PPHQ715_Intercept,PPRQ1110_Intercept)  1.00     2955     3515
## cor(PPHQ115_Intercept,PPRQ2110_Intercept)  1.00     3168     2662
## cor(PPHQ215_Intercept,PPRQ2110_Intercept)  1.00     3490     3064
## cor(PPHQ315_Intercept,PPRQ2110_Intercept)  1.00     3221     2947
## cor(PPHQ415_Intercept,PPRQ2110_Intercept)  1.00     3707     2714
## cor(PPHQ515_Intercept,PPRQ2110_Intercept)  1.00     3876     2685
## cor(PPHQ615_Intercept,PPRQ2110_Intercept)  1.00     3591     3398
## cor(PPHQ715_Intercept,PPRQ2110_Intercept)  1.00     2878     3219
## cor(PPRQ1110_Intercept,PPRQ2110_Intercept) 1.00     3342     3275
## cor(PPHQ115_Intercept,PPRQ3110_Intercept)  1.00     3294     3003
## cor(PPHQ215_Intercept,PPRQ3110_Intercept)  1.00     3796     3309
## cor(PPHQ315_Intercept,PPRQ3110_Intercept)  1.00     3543     3027
## cor(PPHQ415_Intercept,PPRQ3110_Intercept)  1.00     4066     3592
## cor(PPHQ515_Intercept,PPRQ3110_Intercept)  1.00     3662     2937
## cor(PPHQ615_Intercept,PPRQ3110_Intercept)  1.00     3618     3435
## cor(PPHQ715_Intercept,PPRQ3110_Intercept)  1.00     3018     3192
## cor(PPRQ1110_Intercept,PPRQ3110_Intercept) 1.00     3014     2961
## cor(PPRQ2110_Intercept,PPRQ3110_Intercept) 1.00     3166     3395
## cor(PPHQ115_Intercept,PPOQ117_Intercept)   1.00     4740     3062
## cor(PPHQ215_Intercept,PPOQ117_Intercept)   1.00     4391     3056
## cor(PPHQ315_Intercept,PPOQ117_Intercept)   1.00     3959     3003
## cor(PPHQ415_Intercept,PPOQ117_Intercept)   1.00     3857     3560
## cor(PPHQ515_Intercept,PPOQ117_Intercept)   1.00     3668     3011
## cor(PPHQ615_Intercept,PPOQ117_Intercept)   1.00     4048     3457
## cor(PPHQ715_Intercept,PPOQ117_Intercept)   1.00     3021     3201
## cor(PPRQ1110_Intercept,PPOQ117_Intercept)  1.00     3153     3256
## cor(PPRQ2110_Intercept,PPOQ117_Intercept)  1.00     3133     3488
## cor(PPRQ3110_Intercept,PPOQ117_Intercept)  1.00     3134     3138
## 
## Population-Level Effects: 
##                                       Estimate Est.Error l-95% CI u-95% CI Rhat
## PPHQ115_Intercept[1]                     -6.41      3.83   -13.85     0.87 1.00
## PPHQ115_Intercept[2]                     -2.31      3.69    -9.36     4.91 1.00
## PPHQ215_Intercept[1]                     -1.27      4.07    -9.44     6.48 1.00
## PPHQ215_Intercept[2]                      4.08      3.94    -3.75    11.57 1.00
## PPHQ215_Intercept[3]                      8.87      4.13     0.72    17.12 1.00
## PPHQ315_Intercept[1]                    -10.62      4.16   -19.09    -2.78 1.00
## PPHQ315_Intercept[2]                     -6.35      3.84   -14.26     0.88 1.00
## PPHQ315_Intercept[3]                     -2.29      3.52    -9.32     4.55 1.00
## PPHQ315_Intercept[4]                      0.03      3.65    -7.17     7.19 1.00
## PPHQ415_Intercept[1]                     -5.71      4.10   -14.04     2.00 1.00
## PPHQ415_Intercept[2]                     -2.47      3.91   -10.20     5.06 1.00
## PPHQ415_Intercept[3]                      0.48      3.82    -7.01     7.83 1.00
## PPHQ415_Intercept[4]                      3.78      4.01    -3.73    11.56 1.00
## PPHQ515_Intercept[1]                      0.27      3.73    -7.14     7.47 1.00
## PPHQ515_Intercept[2]                      3.25      3.72    -4.17    10.52 1.00
## PPHQ515_Intercept[3]                      4.66      3.76    -2.85    12.07 1.00
## PPHQ515_Intercept[4]                      8.87      4.12     0.73    17.33 1.00
## PPHQ615_Intercept[1]                     -5.94      3.87   -14.04     1.44 1.00
## PPHQ615_Intercept[2]                     -1.77      3.62    -9.13     5.21 1.00
## PPHQ615_Intercept[3]                      0.87      3.58    -6.23     7.79 1.00
## PPHQ615_Intercept[4]                      4.04      3.82    -3.38    11.63 1.00
## PPHQ715_Intercept[1]                      2.90      3.40    -4.00     9.60 1.00
## PPHQ715_Intercept[2]                      4.69      3.44    -2.02    11.41 1.00
## PPHQ715_Intercept[3]                      6.88      3.52    -0.01    13.82 1.00
## PPHQ715_Intercept[4]                      9.62      3.66     2.58    16.76 1.00
## PPRQ1110_Intercept[1]                    -9.50      3.73   -16.78    -2.00 1.00
## PPRQ1110_Intercept[2]                    -7.00      3.59   -14.07     0.35 1.00
## PPRQ1110_Intercept[3]                    -5.97      3.58   -13.03     1.32 1.00
## PPRQ1110_Intercept[4]                    -4.22      3.56   -11.22     2.93 1.00
## PPRQ1110_Intercept[5]                    -0.99      3.56    -7.95     6.27 1.00
## PPRQ1110_Intercept[6]                     0.49      3.63    -6.48     8.05 1.00
## PPRQ2110_Intercept[1]                   -16.36      4.53   -25.47    -7.86 1.00
## PPRQ2110_Intercept[2]                   -12.87      4.23   -21.45    -4.67 1.00
## PPRQ2110_Intercept[3]                   -11.09      4.05   -19.32    -3.11 1.00
## PPRQ2110_Intercept[4]                   -10.18      4.01   -18.19    -2.26 1.00
## PPRQ2110_Intercept[5]                    -6.47      3.85   -14.23     1.31 1.00
## PPRQ2110_Intercept[6]                    -3.56      3.95   -11.43     4.31 1.00
## PPRQ3110_Intercept[1]                   -14.54      4.06   -22.54    -6.66 1.00
## PPRQ3110_Intercept[2]                   -11.35      3.78   -18.87    -4.01 1.00
## PPRQ3110_Intercept[3]                    -9.70      3.61   -16.75    -2.53 1.00
## PPRQ3110_Intercept[4]                    -8.93      3.58   -15.93    -1.82 1.00
## PPRQ3110_Intercept[5]                    -6.78      3.51   -13.59     0.34 1.00
## PPRQ3110_Intercept[6]                    -3.55      3.58   -10.31     3.82 1.00
## PPOQ117_Intercept[1]                     -8.64      3.70   -16.12    -1.60 1.00
## PPOQ117_Intercept[2]                     -7.63      3.60   -14.67    -0.77 1.00
## PPOQ117_Intercept[3]                     -4.34      3.36   -10.81     2.33 1.00
## PPOQ117_Intercept[4]                     -2.37      3.38    -8.87     4.23 1.00
## PPHQ115_Age                               1.35      0.96    -0.53     3.22 1.00
## PPHQ115_Gender01                         -1.23      1.36    -3.92     1.46 1.00
## PPHQ115_Occupation                       -1.87      2.01    -5.84     2.17 1.00
## PPHQ115_Living14Iliveinsharedhousing      2.38      1.90    -1.23     6.27 1.00
## PPHQ115_Living14Ilivewithapartner         0.87      1.52    -2.20     3.93 1.00
## PPHQ115_Living14Ilivewithmyfamily         0.23      1.52    -2.82     3.22 1.00
## PPHQ115_t                                -0.61      0.92    -2.40     1.12 1.00
## PPHQ215_Age                              -3.71      1.15    -6.08    -1.51 1.00
## PPHQ215_Gender01                          1.26      1.50    -1.70     4.30 1.00
## PPHQ215_Occupation                        2.26      2.13    -1.92     6.40 1.00
## PPHQ215_Living14Iliveinsharedhousing     -0.45      1.92    -4.06     3.33 1.00
## PPHQ215_Living14Ilivewithapartner        -0.92      1.67    -4.26     2.36 1.00
## PPHQ215_Living14Ilivewithmyfamily         1.45      1.63    -1.74     4.61 1.00
## PPHQ215_t                                -2.01      0.99    -4.00    -0.13 1.00
## PPHQ315_Age                              -2.17      0.95    -4.05    -0.31 1.00
## PPHQ315_Gender01                         -5.48      1.55    -8.67    -2.55 1.00
## PPHQ315_Occupation                        2.28      2.00    -1.77     6.18 1.00
## PPHQ315_Living14Iliveinsharedhousing     -1.18      1.72    -4.71     2.19 1.00
## PPHQ315_Living14Ilivewithapartner         2.17      1.55    -0.71     5.30 1.00
## PPHQ315_Living14Ilivewithmyfamily        -0.80      1.41    -3.59     1.94 1.00
## PPHQ315_t                                -1.24      0.88    -3.00     0.46 1.00
## PPHQ415_Age                              -0.85      0.89    -2.61     0.95 1.00
## PPHQ415_Gender01                         -3.22      1.45    -6.10    -0.46 1.00
## PPHQ415_Occupation                        3.09      2.09    -0.92     7.14 1.00
## PPHQ415_Living14Iliveinsharedhousing     -1.36      1.89    -5.19     2.33 1.00
## PPHQ415_Living14Ilivewithapartner        -0.97      1.51    -4.04     2.01 1.00
## PPHQ415_Living14Ilivewithmyfamily        -1.41      1.46    -4.16     1.63 1.00
## PPHQ415_t                                -0.66      0.83    -2.33     1.00 1.00
## PPHQ515_Age                              -1.77      0.86    -3.46    -0.08 1.00
## PPHQ515_Gender01                          0.49      1.34    -2.11     3.12 1.00
## PPHQ515_Occupation                        2.31      2.03    -1.76     6.36 1.00
## PPHQ515_Living14Iliveinsharedhousing     -2.72      1.82    -6.32     0.63 1.00
## PPHQ515_Living14Ilivewithapartner        -1.94      1.43    -4.72     0.88 1.00
## PPHQ515_Living14Ilivewithmyfamily        -1.92      1.41    -4.70     0.93 1.00
## PPHQ515_t                                 0.01      0.83    -1.60     1.63 1.00
## PPHQ615_Age                              -1.12      0.85    -2.81     0.54 1.00
## PPHQ615_Gender01                         -2.51      1.31    -5.14     0.02 1.00
## PPHQ615_Occupation                        2.88      2.02    -1.12     6.78 1.00
## PPHQ615_Living14Iliveinsharedhousing     -1.56      1.79    -5.21     1.88 1.00
## PPHQ615_Living14Ilivewithapartner        -0.96      1.47    -3.84     1.98 1.00
## PPHQ615_Living14Ilivewithmyfamily        -0.63      1.41    -3.51     2.04 1.00
## PPHQ615_t                                -1.71      0.89    -3.47    -0.02 1.00
## PPHQ715_Age                              -1.35      0.75    -2.86     0.14 1.00
## PPHQ715_Gender01                         -0.63      1.28    -3.11     1.94 1.00
## PPHQ715_Occupation                        3.80      1.89     0.14     7.49 1.00
## PPHQ715_Living14Iliveinsharedhousing     -2.18      1.75    -5.65     1.14 1.00
## PPHQ715_Living14Ilivewithapartner        -0.39      1.33    -2.97     2.28 1.00
## PPHQ715_Living14Ilivewithmyfamily        -3.08      1.32    -5.65    -0.53 1.00
## PPHQ715_t                                 0.05      0.85    -1.63     1.79 1.00
## PPRQ1110_Age                             -0.17      0.74    -1.66     1.24 1.00
## PPRQ1110_Gender01                         0.43      1.20    -1.95     2.82 1.00
## PPRQ1110_Occupation                      -3.15      1.93    -6.85     0.69 1.00
## PPRQ1110_Living14Iliveinsharedhousing    -2.62      1.87    -6.21     1.14 1.00
## PPRQ1110_Living14Ilivewithapartner        2.58      1.37    -0.06     5.32 1.00
## PPRQ1110_Living14Ilivewithmyfamily       -0.24      1.20    -2.69     2.10 1.00
## PPRQ1110_t                               -0.02      0.77    -1.56     1.50 1.00
## PPRQ2110_Age                              1.75      0.93    -0.06     3.62 1.00
## PPRQ2110_Gender01                        -2.75      1.34    -5.46    -0.20 1.00
## PPRQ2110_Occupation                      -2.45      2.03    -6.57     1.56 1.00
## PPRQ2110_Living14Iliveinsharedhousing     0.37      1.93    -3.33     4.25 1.00
## PPRQ2110_Living14Ilivewithapartner        1.59      1.52    -1.46     4.51 1.00
## PPRQ2110_Living14Ilivewithmyfamily       -2.70      1.50    -5.74     0.10 1.00
## PPRQ2110_t                               -0.43      0.82    -2.04     1.12 1.00
## PPRQ3110_Age                              2.71      0.86     0.94     4.37 1.00
## PPRQ3110_Gender01                        -1.23      1.25    -3.70     1.21 1.00
## PPRQ3110_Occupation                      -3.85      1.95    -7.67     0.05 1.00
## PPRQ3110_Living14Iliveinsharedhousing     0.49      1.75    -2.92     4.01 1.00
## PPRQ3110_Living14Ilivewithapartner        1.32      1.40    -1.43     4.06 1.00
## PPRQ3110_Living14Ilivewithmyfamily       -2.00      1.36    -4.75     0.58 1.00
## PPRQ3110_t                               -1.19      0.82    -2.78     0.38 1.00
## PPOQ117_Age                               2.16      0.87     0.47     3.90 1.00
## PPOQ117_Gender01                         -0.08      1.17    -2.38     2.26 1.00
## PPOQ117_Occupation                       -2.46      1.88    -6.21     1.26 1.00
## PPOQ117_Living14Iliveinsharedhousing      1.00      1.68    -2.27     4.34 1.00
## PPOQ117_Living14Ilivewithapartner         1.41      1.43    -1.39     4.24 1.00
## PPOQ117_Living14Ilivewithmyfamily        -1.22      1.29    -3.80     1.28 1.00
## PPOQ117_t                                -1.68      0.84    -3.39    -0.08 1.00
##                                       Bulk_ESS Tail_ESS
## PPHQ115_Intercept[1]                      3603     2972
## PPHQ115_Intercept[2]                      3524     3048
## PPHQ215_Intercept[1]                      3490     3380
## PPHQ215_Intercept[2]                      3258     3303
## PPHQ215_Intercept[3]                      3243     3015
## PPHQ315_Intercept[1]                      3034     3174
## PPHQ315_Intercept[2]                      3570     3301
## PPHQ315_Intercept[3]                      3673     3378
## PPHQ315_Intercept[4]                      3747     2978
## PPHQ415_Intercept[1]                      3052     3072
## PPHQ415_Intercept[2]                      3266     2915
## PPHQ415_Intercept[3]                      3405     2802
## PPHQ415_Intercept[4]                      3344     2788
## PPHQ515_Intercept[1]                      3655     3265
## PPHQ515_Intercept[2]                      3718     3268
## PPHQ515_Intercept[3]                      3722     3454
## PPHQ515_Intercept[4]                      3560     3214
## PPHQ615_Intercept[1]                      3507     3217
## PPHQ615_Intercept[2]                      3967     3342
## PPHQ615_Intercept[3]                      4194     3333
## PPHQ615_Intercept[4]                      4269     3467
## PPHQ715_Intercept[1]                      3600     2977
## PPHQ715_Intercept[2]                      3524     2962
## PPHQ715_Intercept[3]                      3435     3143
## PPHQ715_Intercept[4]                      3210     3239
## PPRQ1110_Intercept[1]                     3729     3222
## PPRQ1110_Intercept[2]                     3982     3322
## PPRQ1110_Intercept[3]                     3896     3126
## PPRQ1110_Intercept[4]                     3938     3431
## PPRQ1110_Intercept[5]                     3745     3454
## PPRQ1110_Intercept[6]                     3611     3426
## PPRQ2110_Intercept[1]                     3017     2799
## PPRQ2110_Intercept[2]                     3607     2577
## PPRQ2110_Intercept[3]                     3812     2457
## PPRQ2110_Intercept[4]                     3814     2719
## PPRQ2110_Intercept[5]                     3768     2661
## PPRQ2110_Intercept[6]                     3875     2814
## PPRQ3110_Intercept[1]                     2717     2729
## PPRQ3110_Intercept[2]                     2983     2986
## PPRQ3110_Intercept[3]                     2977     3054
## PPRQ3110_Intercept[4]                     2958     3066
## PPRQ3110_Intercept[5]                     2969     2775
## PPRQ3110_Intercept[6]                     3125     3215
## PPOQ117_Intercept[1]                      3224     3114
## PPOQ117_Intercept[2]                      3363     3164
## PPOQ117_Intercept[3]                      3762     3264
## PPOQ117_Intercept[4]                      3779     3205
## PPHQ115_Age                               3724     2640
## PPHQ115_Gender01                          3579     3424
## PPHQ115_Occupation                        3197     3160
## PPHQ115_Living14Iliveinsharedhousing      3866     2994
## PPHQ115_Living14Ilivewithapartner         3250     3081
## PPHQ115_Living14Ilivewithmyfamily         3838     3302
## PPHQ115_t                                 7583     3130
## PPHQ215_Age                               2644     2832
## PPHQ215_Gender01                          2724     2652
## PPHQ215_Occupation                        3183     3077
## PPHQ215_Living14Iliveinsharedhousing      3522     3005
## PPHQ215_Living14Ilivewithapartner         3121     2933
## PPHQ215_Living14Ilivewithmyfamily         3267     3285
## PPHQ215_t                                 4685     2920
## PPHQ315_Age                               2648     2837
## PPHQ315_Gender01                          3015     3430
## PPHQ315_Occupation                        3260     3000
## PPHQ315_Living14Iliveinsharedhousing      4254     2844
## PPHQ315_Living14Ilivewithapartner         2926     3380
## PPHQ315_Living14Ilivewithmyfamily         3636     3253
## PPHQ315_t                                 6341     2827
## PPHQ415_Age                               3043     2726
## PPHQ415_Gender01                          3843     2840
## PPHQ415_Occupation                        3324     2937
## PPHQ415_Living14Iliveinsharedhousing      3755     2774
## PPHQ415_Living14Ilivewithapartner         3104     3044
## PPHQ415_Living14Ilivewithmyfamily         2774     3288
## PPHQ415_t                                 6508     2787
## PPHQ515_Age                               3181     2669
## PPHQ515_Gender01                          3954     3543
## PPHQ515_Occupation                        3140     2918
## PPHQ515_Living14Iliveinsharedhousing      4105     3313
## PPHQ515_Living14Ilivewithapartner         3234     3032
## PPHQ515_Living14Ilivewithmyfamily         3171     2723
## PPHQ515_t                                 6516     3114
## PPHQ615_Age                               3269     3017
## PPHQ615_Gender01                          3877     3363
## PPHQ615_Occupation                        3392     2693
## PPHQ615_Living14Iliveinsharedhousing      4413     3509
## PPHQ615_Living14Ilivewithapartner         3833     2873
## PPHQ615_Living14Ilivewithmyfamily         3569     2914
## PPHQ615_t                                 5729     3198
## PPHQ715_Age                               2977     2737
## PPHQ715_Gender01                          4099     2948
## PPHQ715_Occupation                        2817     2731
## PPHQ715_Living14Iliveinsharedhousing      4163     3296
## PPHQ715_Living14Ilivewithapartner         4276     2841
## PPHQ715_Living14Ilivewithmyfamily         3578     3012
## PPHQ715_t                                 6925     2766
## PPRQ1110_Age                              3211     3171
## PPRQ1110_Gender01                         4291     3063
## PPRQ1110_Occupation                       3386     3140
## PPRQ1110_Living14Iliveinsharedhousing     4294     3427
## PPRQ1110_Living14Ilivewithapartner        3510     2865
## PPRQ1110_Living14Ilivewithmyfamily        3758     3193
## PPRQ1110_t                                6460     3176
## PPRQ2110_Age                              3311     3330
## PPRQ2110_Gender01                         3526     3295
## PPRQ2110_Occupation                       3531     2934
## PPRQ2110_Living14Iliveinsharedhousing     4357     3126
## PPRQ2110_Living14Ilivewithapartner        3349     2790
## PPRQ2110_Living14Ilivewithmyfamily        2901     2965
## PPRQ2110_t                                7172     2682
## PPRQ3110_Age                              2825     3080
## PPRQ3110_Gender01                         3806     3303
## PPRQ3110_Occupation                       2818     2956
## PPRQ3110_Living14Iliveinsharedhousing     4061     3012
## PPRQ3110_Living14Ilivewithapartner        2950     2948
## PPRQ3110_Living14Ilivewithmyfamily        3365     2714
## PPRQ3110_t                                5891     3233
## PPOQ117_Age                               3038     3084
## PPOQ117_Gender01                          3441     3217
## PPOQ117_Occupation                        3281     2729
## PPOQ117_Living14Iliveinsharedhousing      4490     2920
## PPOQ117_Living14Ilivewithapartner         3357     2852
## PPOQ117_Living14Ilivewithmyfamily         3541     2901
## PPOQ117_t                                 4846     3013
## 
## Family Specific Parameters: 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc_PPHQ115      1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPHQ215      1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPHQ315      1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPHQ415      1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPHQ515      1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPHQ615      1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPHQ715      1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPRQ1110     1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPRQ2110     1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPRQ3110     1.00      0.00     1.00     1.00 1.00     4000     4000
## disc_PPOQ117      1.00      0.00     1.00     1.00 1.00     4000     4000
## 
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Gender has a negative impact (male/transmale answers lower on the response), however it is not significant.

Living with a partner makes respondents answer with higher response in Q11, however it is not significant.

4 Computational environment

print(sessionInfo(), locale=FALSE)
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] latex2exp_0.4.0   ggthemes_4.2.4    here_1.0.1        data.table_1.14.0
##  [5] openxlsx_4.2.3    patchwork_1.1.1   bayesplot_1.8.0   forcats_0.5.1    
##  [9] stringr_1.4.0     dplyr_1.0.4       purrr_0.3.4       readr_1.4.0      
## [13] tidyr_1.1.2       tibble_3.1.0      ggplot2_3.3.3     tidyverse_1.3.0  
## [17] plyr_1.8.6        brms_2.14.11      Rcpp_1.0.6       
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1            backports_1.2.1         igraph_1.2.6           
##   [4] splines_4.0.4           crosstalk_1.1.1         TH.data_1.0-10         
##   [7] rstantools_2.1.1        inline_0.3.17           digest_0.6.27          
##  [10] htmltools_0.5.1.1       rsconnect_0.8.16        fansi_0.4.2            
##  [13] magrittr_2.0.1          modelr_0.1.8            RcppParallel_5.0.3     
##  [16] matrixStats_0.58.0      xts_0.12.1              sandwich_3.0-0         
##  [19] prettyunits_1.1.1       colorspace_2.0-0        rvest_0.3.6            
##  [22] haven_2.3.1             xfun_0.21               callr_3.5.1            
##  [25] crayon_1.4.1            jsonlite_1.7.2          lme4_1.1-26            
##  [28] survival_3.2-7          zoo_1.8-8               glue_1.4.2             
##  [31] gtable_0.3.0            emmeans_1.5.4           V8_3.4.0               
##  [34] pkgbuild_1.2.0          rstan_2.26.0.9000       abind_1.4-5            
##  [37] scales_1.1.1            mvtnorm_1.1-1           DBI_1.1.1              
##  [40] miniUI_0.1.1.1          xtable_1.8-4            stats4_4.0.4           
##  [43] StanHeaders_2.26.0.9000 DT_0.17                 htmlwidgets_1.5.3      
##  [46] httr_1.4.2              threejs_0.3.3           ellipsis_0.3.1         
##  [49] pkgconfig_2.0.3         loo_2.4.1               sass_0.3.1             
##  [52] dbplyr_2.1.0            utf8_1.1.4              tidyselect_1.1.0       
##  [55] rlang_0.4.10            reshape2_1.4.4          later_1.1.0.1          
##  [58] munsell_0.5.0           cellranger_1.1.0        tools_4.0.4            
##  [61] cli_2.3.1               generics_0.1.0          broom_0.7.5            
##  [64] ggridges_0.5.3          evaluate_0.14           fastmap_1.1.0          
##  [67] yaml_2.2.1              processx_3.4.5          knitr_1.31             
##  [70] fs_1.5.0                zip_2.1.1               nlme_3.1-152           
##  [73] mime_0.10               projpred_2.0.2          xml2_1.3.2             
##  [76] rstudioapi_0.13         compiler_4.0.4          shinythemes_1.2.0      
##  [79] curl_4.3                gamm4_0.2-6             reprex_1.0.0           
##  [82] statmod_1.4.35          bslib_0.2.4             stringi_1.5.3          
##  [85] ps_1.5.0                Brobdingnag_1.2-6       lattice_0.20-41        
##  [88] Matrix_1.3-2            nloptr_1.2.2.2          markdown_1.1           
##  [91] shinyjs_2.0.0           vctrs_0.3.6             pillar_1.5.0           
##  [94] lifecycle_1.0.0         jquerylib_0.1.3         bridgesampling_1.0-0   
##  [97] estimability_1.3        httpuv_1.5.5            R6_2.5.0               
## [100] bookdown_0.21           promises_1.2.0.1        gridExtra_2.3          
## [103] codetools_0.2-18        boot_1.3-27             colourpicker_1.1.0     
## [106] MASS_7.3-53.1           gtools_3.8.2            assertthat_0.2.1       
## [109] rprojroot_2.0.2         withr_2.4.1             shinystan_2.5.0        
## [112] multcomp_1.4-16         mgcv_1.8-34             parallel_4.0.4         
## [115] hms_1.0.0               grid_4.0.4              coda_0.19-4            
## [118] minqa_1.2.4             rmarkdown_2.7           cmdstanr_0.3.0.9000    
## [121] shiny_1.6.0             lubridate_1.7.9.2       base64enc_0.1-3        
## [124] dygraphs_1.1.1.6

5 References